Macroscopic quantum bound states of Bose Einstein condensates in optical lattices 
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We discuss localized ground states of the periodic Gross-Pitaevskii equation in the framework of a 
quantum linear Schrodinger equation with effective potential determined in self-consistent manner. 
We show that depending on the interaction among the atoms being attractive or repulsive, bound 
states of the linear self consistent problem are formed in the forbidden zones of the linear spectrum 
below or above the energy bands. These eigenstates are shown to be exact solitons of the GPE 
equation. The implications of this bound state interpretation on the existence of a delocalization 
transition for multidimensional solitons is briefly discussed. 
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One interesting phenomenon occurring in periodic non- 
linear systems is the possibility to stabilize localized ex- 
citations as a result of the interplay between periodicity 
and nonlinearity. An example of this is provided by the 
nonlinear Schrodinger equation (NLS) with periodic po- 
tential. It is well known that the defocusing NLS does 
not admit bright soliton solutions, these being unstable 
against background decay [1]. The presence of a peri- 
odic potential, however, allows to stabilize bright soli- 
tons against decay, a phenomenon which is presently in- 
vestigated in connection with Bose Einstein condensates 
(BEG) in optical lattices (OL). The possibility to form 
bright solitons in repulsive BEG with OL was analyti- 
cally and numerically demonstrated, both for a discrete 
version of the NLS describing BEG arrays in the tight- 
binding approximation [2] and for the Gross-Pitaevskii 
equation (GPE) describing the properties of a continuous 
BEG in the mean field approximation [3-5] . The mecha- 
nism underlying soliton formation in periodic structures 
was identified to be the modulational instability of the 
Bloch states at the edges of the Brillouin zone [4] . These 
localized excitations correspond to states with energies 
inside the gaps of the underlying linear band structure 
(in nonlinear optics they are called gap solitons) and with 
an effective mass which depends on the sign of the in- 
teraction (for repulsive interactions, bright solitons have 
negative effective mass, this explaining their existence in 
BEG with OL [4,6]). The usage of hnear concepts such 
as Bloch states, effective mass, etc. [4,6,7], makes natural 
to ask whether nonlinear states could be interpreted in a 
pure linear (quantum mechanical) context. 

The aim of the present paper is to address this prob- 
lem by showing that soliton solutions of the periodic non- 
linear Schrodinger equations correspond to bound states 
of the linear Schrodinger equation with an effective po- 
tential which can be determined in self-consistent (SG) 
manner. This problem will be discussed on the physical 
example of a Bose Einstein condensate in an optical lat- 
tice (OL) described, in mean field approximation, by the 
following normalized Gross-Pitaevskii equation 



where x is the nonlinear parameter, x denotes three di- 
mensional coordinates and J7(x) is a periodic potential 
representing the OL. To discuss bound state features of 
solitons we restrict to the one dimensional case (the ap- 
proach however is of general validity and can be applied 
to NLS type equations in arbitrary dimensions). At the 
end of the paper we will briefly discuss the implications of 
the bound state interpretation of localized solutions on 
the soliton delocalization transition observed in higher 
dimensions [8]. We remark that the properties of soli- 
tons of the GPE in optical lattices were studied in [5] 
in terms of orbits of a chaotic system. Self-consistent 
approaches were also used as numerical tools to study 
discrete breathers of the discrete NLS [10] and the sta- 
bility of gap solitons [9] . The physical implications and 
the full potentiality of the SG approach, however, have 
not been investigated. 

Our analysis is based on the simple observation 
that the stationary localized ground states ips{x,t) — 
^{x) exp(— /ii) of the GPE (and more generally of any 
nonlinear Schrodinger-like equation) can be obtained by 
solving in a self-consistent manner the following linear 
Schrodinger problem 



(2) 



(1) 



with the effective potential 

Veff = Uoi{x) + Us{x) = Acos{2x)+x\i^s{x)\^- (3) 

Here Uoi = v4cos(2a;) is the OL and Us is the potential 
associated with a given eigenstate of the quantum prob- 
lem (2). For a self-consistent solution, one starts with 
a trial wavefunction for V's (typically a gaussian wave- 
form), calculates the effective potential and solves the 
corresponding eigenvalue problem (2). Then, one selects 
a given eigenstate (for example the ground state but not 
necessarily) as new trial function and iterates the proce- 
dure until convergence is reached. 
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FIG. 1. Panel (a) Energy spectrum for the effective po- 
tential (3) with A = 3 and x = (Mathieu equation). Full 
lines represent exact values of the band edges of the Math- 
ieu equation while dots arc the eigenvalues obtained with the 
above procedure on a lattice of length L = 40n, with A'^ = 512 
points. Panel (b) Lowest energy band for the effective poten- 
tial in Eq. (3) with tps taken as the ground state of the system 
and for X ~ ~^ (attractive case). Parameters are fixed as in 
panel (a). Panel (c) The same as in panel (b) but for A = —3. 
Panel (d) Transition from the metastable IS mode to the OS 
ground state corresponding to the lower level of panel (c). The 
optical lattice (scaled by a factor 3) is reported as an help to 
locate the symmetry center of the solutions. Parameters are 
fixed as in panel (c). 

The problem is thus reduced to the diagonalization of 
the quantum Hamiltonian 



H = Ho + Veff{x) 



(4) 



with Ho = — the kinetic energy operator. This can be 
effectively done by adopting a discrete coordinate space 
representation = 77.0,}, n = l,...,N, with a = L/N 
the discretization constant, L the size of the system and 
N the total number of points. A basis for this space is 



simply a basis of R 



N 



the set of N-component vec- 



tors of the type \n) = (0, ...0, 1, 0, 0), with the 1 in the 
position n. The effective potential is obviously diagonal 
in this representation i.e. (77|V^;//|r7') — Vef f{na)dn.n' , 
while Hq is diagonal in the reciprocal representation \kn), 
(kn = 211 /Ln, the two representations being related by 
the Fourier transform (unitary transformation). The ma- 
trix elements of the Hamiltonian H can then be con- 
structed as 

{n\H\n') = Hn,n' = {n\F-^ HoF\n') + V,ff{na)Sn,n' (5) 



where F\n) denotes the Fourier transform of the vector 
In). 
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FIG. 2. Wavefunctions and corresponding effective poten- 
tial for the bound states below the lowest bands of Fig la for 
attractive interaction ^ = — 1- Panel (a). OS mode and cor- 
responding effective potential for E — —1.1667950 (ground 
state) and A = 3. The effective potential was scaled by a 
factor 6 for graphical convenience. Panel (b). Same as (a) 
for the IS mode at £; = -1.0485745 and A = -3. Panel (c). 
Same as (a) for the OA mode. E = -0.999261. Panel (d). 
Same as (b) for the lA mode. E = -0.9947127. (e) Energy 
levels of the OS, IS, OA, I A, nodes inside the gap between the 
first two bands. The continuous line denotes the lower edge of 
the second band of Fig. la while the dotted lines denote the 
degenerate levels. Parameters are fixed eis for corresponding 
modes in panels a-d. (f). Wavefunctions associated to the 
levels in panel e. For graphical convenience the IS mode Wcis 
shifted by -1.0 down while the I A and OS modes were shifted 
up by 0.5 and 1.0, respectively. 

For an effective construction of these matrix elements 
one can use the fast Fourier transform while for the 
computation of the spectrum one can recourse to stan- 
dard numerical routines for the diagonalization of real 
symmetric matrices. To check the method wc consider 
first the case of a linear effective potential of the form 
Veff — Uoi — Acos(2a;) for which the eigenvalue prob- 
lem reduces to the wcill known Mathieu equation. In Fig. 
la we depict the lowest part of the spectrum (notice that 
in this case there is no SC procedure due to the linear- 
ity of the problem) from which we see the appearance 
of a band structure with band edges which exactly coin- 
cide with the values obtained for the Mathieu equation 
(for high energy bands to get good accuracy one needs 
to increase N). In this paper we are mainly interested in 
the localized states associated with the lowest two bands 
(i.e., the ones physically most relevant), and for this pur- 
pose the choice of N = 256 will be adequate for most of 
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the calculations. 



decay). 
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FIG. 3. Panel (a). Time evolution of the OS bound state 
of Fig 2a as obtained from GPE. Panel (b). Same as in panel 
(a) for the IS mode of Fig. 2c. 
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FIG. 4. Energy spectrum in correspondence of the local- 
ized states above the lowest band of Fig. la for repulsive 
interactions x = 1- Panel (a). Spectrum associated to the 
OS mode. Parameters arc A = —3, A'' = 256, L = 40-k. Panel 
(b). Same as panel (a) but for the lA mode with A = 3. 
The continuous lines denote exact band edges of the Mathieu 
equation. 

In Fig. lb we show how the lowest band of panel 

la is modified by the nonlinear potential Veff{x) = 
Acos(2x) + xlV'oP, where ipo is taken to be the ground 
state of the system, for the case x < (negative scatter- 
ing length). A bound state below the band which rapidly 
converges to a constant value is quite evident. One can 
see that the state forms from the lower edge of the band 
and is accompanied by a rearrangement of the extended 
states inside the band. The corresponding eigenvector 
is depicted in Fig. 2a together with its effective potential. 
Notice that the potential has an attractive character (po- 
tential well) and the bound state is symmetric around a 
minima of the OL, i.e., it resembles the onsite-symmetric 
intrinsic localized mode (ILM) of nonlinear lattices (NL) 
[11]. In the following we shall call it the onsite symmet- 
ric (OS) bound state. By shifting the phase of the OL 
by TT (i.e. by changing the sign of A) one obtains an 
eigenstate centered on maxima instead than on minima. 
This bound state is depicted in Fig. 2b and in anal- 
ogy with NL we shall call it the intersite symmetric (IS) 
mode. The corresponding spectrum is reported in Fig. 
Ic. Notice that the IS mode corresponds to the plateau 
formed just before the decay into the OS mode occurs 
as shown in panel Id (also note the appearance of an 
intersite-symmetric (lA) excited level in Fig. Ic which is 
absorbed into the band in correspondence of the IS-OS 
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FIG. 5. Wavefunctions and effective potentials of the 
bound states levels of Fig. 4 a,b, for the repulsive case x = 1- 
Panel (a). OS mode with corresponding effective potential 
(thin line). Energy is E = -0.078355 and A = -3. Panel 
(b). Same as in Panel (a) for the lA mode. E = -0.376645, 
A = 3. Panel (c). Same as in panel (a) for the OA mode. 
E = -0.683070. Panel (d). Same in panel (b) for the IS 
mode. E — —0.691676. Parameters arc fixed as A'" = 256, 
L = 407r for panels (a), (b), and N = 512, L = 407r for pan- 
els (c), (d). The effective potentials have been reduced by a 
factor 6 for graphical convenience. 

We have checked that these bound states coincide with 

the ones obtained with the approach of Ref. [5] for the 
same values of energy. The stability of the OS mode and 
the decay of the IS mode into the OS state was checked 
by direct numerical integrations of the GPE (see Fig. 3). 
To obtain the onsite asymmetric (OA) mode one needs 
to take the first excited state ipi as effective potential in 
the SC procedure. This indeed produces an exact soliton 
solution of the GPE of type OA as shown in Fig. 2c. A 
shifting of the potential by tt produce the intersite asym- 
metric (lA) mode of Fig. 2d. These solutions have the 
same energies and are more unstable than the IS mode 
(they, however, do not decay into the ground state but 
get mixed with the extended states in the band). In gen- 
eral, the effective potentials can be taken of the form 
^e// = ^oi + xlV'n(ic)| with tjjn the n-th eigenstate of the 
eigenvalue problem in (2). If the energy of ipn lies outside 
the bands a localized mode of the type described above 
is produced, while if it lies inside a band, extended states 
which are nonlinear analogue of the Bloch states [7] , are 
produced. From this we conclude that both localized 
and extended solutions of the GPE are exact quantum 
eigenstates of the Schrodinger equation with a suitable 
effective potential. From the above analysis one expects 
that below each higher energy bands, eigenstates of the 
same symmetry type as the ones found for the lowest 
band should exist. This is precisely what we show in 
Fig.s 2e, 2f for the energy spectrum and the correspond- 
ing eigenstated in the gap below the second band. 
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FIG. 6. Panel (a). Two soliton bound state of the repul- 
sive GPE obtained from the SC procedure by using as effective 
potential Veff = Vol + l^i'i — where '4>i,ip2 denote two 
eigenfunctions at the top of the first band. Panel (b). Time 
evolution (modulo square) of the bound state in panel (a) 
taken as initial condition for the integration of the full GPE. 

Notice that the OA and the IS eigenstates are degen- 
erated (the same occurs also to the OS and lA modes). 
The OA and IS bound states are both very stable while 
the energy levels of the OA and IS modes, after estab- 
lishing a plateau similar to the one in Fig Ic, become 
unstable (the energy oscillates between this level and the 
lower edge of the second band). The instability of these 
modes can be understood as an hybridization of the state 
(being very close to the band edge) with extended states 
of the second band and is confirmed by direct numerical 
integration of the PDE system. 

Similar localized modes can be found also for repulsive 
interactions (x > 0), the main diflference with the pre- 
vious case being that now the states appear in the gap 
above the band edges instead than below. This is shown 
in Fig. 4 for the lowest energy levels inside the first gap. 
The corresponding eigenvectors are shown in Fig. 5 to- 
gether with their effective potentials. Notice that the 
potential has a local repulsive character (it increases in 
correspondence of the states) so that these bound states 
could not exist without the OL. We remark that localized 
modes similar to the ones described in this paper were 
found also in atomic-molecular BEG using an approach 
based on Wannier functions [12]. 

Besides localized and extended states, the SC proce- 
dure allows to construct full nonlinear bands in reciprocal 
space (we omit details for brevity). It is worth remarking 
that more complicated set of solutions of the GPE can 
be constructed with the SC procedure by taking as ef- 
fective potentials linear combinations of eigenstates. An 
example of this is shown in Fig. 6 for the case of repul- 
sive interaction. We see that a linear combination of two 
eigenstates leads to a bound state with two humps which 
corresponds to a multi-soliton solution of the GPE (see 
panel b). This is a general property and we conjecture 
that all solutions of the periodic GPE ( and more in gen- 
eral of the NLS-like equations with arbitrary potentials) 



can he obtained with the SC method taking all possible 
combinations of linear eigenstates as effective potentials. 

Before closing this paper we wish to remark that the 
above bound state interpretation has important conse- 
quences on the delocahzing transition [8] of localized so- 
lutions of the GPE in OL. Since in a ID potential bound 
state always exists, it follows from the above analysis that 
no delocahzing transition of a BEG soliton can occur in 
this case. On the contrary, for _D > 2 a finite potential 
depth is required to form a bound state, this implying 
that a soliton delocalization transition can occur. A de- 
tailed investigation of the delocahzing transition of BEG 
solitons in OL will be reported elsewhere [13]. 
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